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Abstract 

We consider a quadrilateral ’mini’ finite element for approximating the solution 
of Stokes equations using a quadrilateral mesh. We use the standard bilinear finite 
element space enriched with element-wise defined bubble functions for the velocity 
and the standard bilinear finite element space for the pressure space. With a simple 
modification of the standard bubble function we show that a single bubble function is 
sufficient to ensure the inf-sup condition. We have thus improved an earlier result on 
the quadrilateral ’mini’ element, where more than one bubble function are used to get 
the stability. 

Index terms — Stokes equations, mixed finite elements, Mini finite element, inf-sup 
condition, bubble function 

AMS subject classification. 65N30, 65N15, 74B10 


1 Introduction 

A very simple finite element method for the Stokes problem for a simplicial mesh is presented 
by Arnold, Brezzi and Frotin pQ, where the velocity space is discretised by using the standard 
linear finite element space enriched with element-wise bubble functions and the pressure 
space is discretised by using the standard linear finite element space. The enrichment of the 
velocity space is done to ensure the stability of the finite element method, and this increases 
one vector degree of freedom per element. An extension of the finite element method to 
the quadrilateral mesh is done by Bai [2] , where the author enriches the velocity space with 
more than a single vector bubble function per element. The inf-sup condition is proved by 
using a macro element technique nsi, where a single element is used as a macro element. 

In this article we show that with a small modification of the standard bubble function 
we can get the stability just by using a single vector bubble function per element. The 
main difference with the technique proposed by Bai j2] is that it is not possible to show 
the inf-sup condition using a single element as a macro element. We need to use a macro 
element consisting of four elements to prove the inf-sup condition in our situation. Another 
relevant finite element method is presented by Lamichhane [8j, where two different meshes 
are used to discretise the velocity and the pressure space, and a single vector bubble degree 
of freedom per element is used to get the stability. The pressure space is discretised by 
the space of piecewise constant functions on the dual mesh. However, the main difficulty 
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of the technique presented by Lamichhane 0 is that the bubble function is obtained by 
multiplying the standard bubble function by the gradient of a bilinear basis function, and 
hence the bubble function cannot be defined on a reference element. The standard bubble 
function on the unit square is the lowest degree polynomial which vanishes on the boundary 
of the square. Here we modify the standard bubble function cm to get stability of the 
numerical scheme by using a single vector bubble function per element with a continuous 
pressure approximation. We also investigate two choices of bubble functions, where both of 
them can be defined on a reference element. Since the first mini finite element is introduced 
for simplicial meshes [T] with a single bubble function per element, this new contribution 
gives a unified framework for quadrilaterals and triangles. The idea can easily be extended 
to the three-dimensional case. 


2 Stokes equations 


This section is devoted to the introduction of the boundary value problem of the Stokes equa¬ 
tions. Let flin R 2 , be a bounded domain with polygonal boundary T. For a prescribed body 
force / £ \L 2 {Cl)] 2 , the Stokes equations with homogeneous Dirichlet boundary condition in 
T reads 

—vAu + Vp = f in Cl 
div u = 0 in ft 


( 2 . 1 ) 


with u = 0 on T, where u is the velocity, p is the pressure, and v denotes the viscosity of 
the fluid. 

Here we use standard notations L 2 (Cl), 1L 1 (H) and Hq(CI) for Sobolev spaces, see [HE! 
for details. Let V := [iL^H)] 2 be the vector Sobolev space with inner product (-,-)i,n and 
norm || • Hqn defined in the standard way: ( u,v)i ^ := Vj)i,n, and the norm being 

induced by this inner product. We also define another subspace M of L 2 (Cl) as 


P = £ L 2 (Cl) : J qdx = o| . 


The weak formulation of the Stokes equations is to find (it,p) £ V x P such that 

v f n Vu : Vv dx + div v p dx =£(v), v£V, 
f Q div uqdx = 0, q £ P, 


( 2 . 2 ) 


where £(v) = j n f ■ v dx. It is well-known that the weak formulation of the Stokes problem is 
well-posed [7]. In fact, if the domain Cl is convex, and f £ [L 2 (Cl)] 2 , we have u £ [H 2 (Cl)} 2 , 
p £ H 1 {Cl) and the a priori estimate holds 

IMko + lbl|i,n < C'll/llo.n, 
where the constant C depends on the domain Cl. 


3 Finite element discretizations 

We consider a quasi-uniform triangulation Th of the polygonal domain Cl, where Th consists 
of convex quadrilaterals. The finite element meshes are defined by maps from the reference 
square K = (0, l) 2 . 
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Let Qi(K) be the space of bilinear polynomials in K. We start with the finite element 
space of continuous functions whose restrictions to an element K are obtained by maps of 
bilinear functions from the reference element: 

S h := { v h G v h \ K = v h ° F«\v h G Qi(A'), K G %} , (3.1) 

where Fk : K K is an iso-parametric map. We note that the iso-parametric map Fk is 
generated by using the basis functions of Q\(K). It is clear that if v G Qi(K), then voF^ 1 
is in general not a polynomial on the quadrilateral K. 

In the following we assume that each element K G Th is a parallelogram and the map Fk 
is affine. Let bx be a bi-variate polynomial of x G ffi 2 with bx = 0 on dK and bxixx) = 1 , 
where xx G M 2 is the centroid of K. This is called a bubble function corresponding to the 
element I\ G Th- Defining the space of bubble functions 


Bh := {bh G C7°(n) : bh\x = cxbx , ck G R, K G Th}, 


(3.2) 


we introduce our finite element space for velocity as Vh = [S/j ® B h ] 2 . The finite element 
space for the pressure is taken as the standard bilinear finite element space 

S* h := {rhGLl^nH 1 ^), v h \ K =VhoF^\v h GQ 1 (K), K G %} . (3.3) 


Then, the finite element approximation of (2.21 is defined as a solution to the following 
problem: find ( u h ,p h ) G Vh x S* h such that 


a(u h , v h ) + b(v h ,p h ) = t{v h ), v h G V h , 


Qh e Sv 


(3.4) 


K u h,<lh) = o, 

We need the following conditions to prove that there is a unique solution of the discrete 


problem (3.4) and the discrete solution converges optimally to the continuous solution. 


1. The bilinear forms a(-, •) on Vh x Vh and &(•,•) on Vh x S ^ are continuous. 

2. The bilinear form on Vh x Vh is elliptic. 

3. There exists a constant 0 > 0 independent of the mesh-size such that for any qh G 
we have 


b(v h ,q h ) . on ii 

sup T —tt-> P\\qh\\o,n- 

v h aV h ll^llbn 


(3.5) 


The smallest constant 0 with the property 


0 = inf sup 


b(v h ,qh) 


<JheS* VhG y h Iloilo,n 

is called the inf-sup constant. 


(3.6) 


4 The Macro-Element Technique 


We prove the inf-sup condition (3.5) using a macro-element technique proposed by Stenberg 


m ■ A macro-element M is a connected set of elements in Th- Moreover, two macro-elements 
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M\ and M 2 are said to be equivalent if they can be mapped continuously onto each other 
m- We define the following three spaces associated with the macro-element M,: 

v\ = [H%{Mi)] 2 nV, Si = [v h G H 1 (M i ), v h \ K = v h o F~\v h € Qi(K), I< G T h , K C M,} , 

and 

Bi = {q h G S\\ b(v h , q h ) = 0 , v h G V l h } . 

Moreover, we denote by the set of all edges in 7 ~h interior to The macro-element 
partition A4/, of f l then consists of macro-elements {Mi}fL 1 with Cl = [J^ =1 M,. The macro¬ 
element technique is given by the following theorem 1103- 

Theorem 1 Suppose that there is a fixed set of equivalence classes £j, j = 1 ,■■■ ,q, of 
macro-elements, a positive integer L, and a macro-element partition M h such that 

(Ml) For each Mi G £j, j = 1, • • • ,q, the space Bi is one-dimensional, consisting of func¬ 
tions that are constant on Mi. 


(M2) Each Mi G Mh belongs to one of the classes £j, j = 1, • • • ,q. 

(M3) Each T £ Th is contained in at least one and not more than L macro-elements of A4h- 

(Mf) Each e G Tft, is contained in the interior of at least one and not more than L macro- 
elments of Mi h- 


Then the inf-sup condition (3.51 is satisfied. 



Figure 1: The set Mi, where four elements of Th touch the vertex aq 

In the following we consider a macro-element consisting of four squares as shown in 
Figure [lj With this partition of macro-elements we can see that Assumptions (M2)-(M4) 
are all satisfied. We now show that the proof of Assumption (Ml) depends on the choice of 
bubble functions. 
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4.1 Choice of bubble functions 


For simplicity of calculation we assume that Mj is a parallelogram so that there is an 
invertible affine mapping : S —>• M, ; , which transforms the square S = [—1, l] 2 to Mj with 
the property 


X 

= A t 

£ 

+ 

x 0 

y. 


y. 


yo_ 


(4.1) 


where A t is a 2 by 2 matrix, (x,y) £ Mj and (£, 77 ) £ S. Let VJJ = span{^) fe }^ =1 , V l h = \V£\ 2 
and S l h = span{(jjfc}| =1 . We use the notation fa and (pk to denote corresponding basis 
functions on the square S, where fa and (pk are functions of £ and 77. We have shown the 
numbering of functions fa and ! Pj on the reference square S in Figure pi where we have used 
big circles for the functions in Vh, and small circles for functions in S'l. 

Let Vh £ Vh with Vh = 1 v kfa and Vk £ R 2 . Then 


b(v h ,q h )= / V v h q h dx = y2 v k - \7(f) k q h dx. 

J Mi k=1 J Mi 

Using a chain rule we write 

Vfa=A~ T (vfaoF- 1 ), 

where V denotes the gradient on the reference square S. Let qh = J2j= 1 QjVji an d thus 


I M { 


5 9 

V -v h q h dx = EE qjv k 

k= 1 3 =1 


5 9 


I Mi 


V fa qjipj dx = | clet Ai\ EE QjV k ■ / A i T Vfa<Pj dx. 

k= 1 j =1 


We see that we can find a matrix D such that 



■ v h qh dx 


fDv, 


where 


Qi 


Vi 


Vi 


, and v = 

V2 

= 

V 2 

. 99 . 


V 5 _ 


y w _ 


Thus we need to show that the rank of the matrix D is 8 in order to prove that the dimension 
of the space Bi is one. 

Since Aj is an invertible matrix, the rank of the matrix will be unchanged if we replace Mi 
by the reference element S, so that we want to investigate the rank of the matrix D £ R 10x9 , 
where the jth row of D is 


/ d^faipj dx, / d n fa(pj dx, / (\faPj dx, / d v fa fa dx, ■ ■ ■ , / dx, / d v fa fa dx 
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Figure 2: The numbering of functions <j>k and (f>j on the reference square S 

4.1.1 Standard bubble functions 

Consider the unit square K = (0, l) 2 in two dimensions. We start with the standard choice 
of the bubble function bj< = \0>xy(l — x)(l — y). The matrix D is explicitly computed as 



We compute the rank of this matrix using MAPLE and obtain it to be 7. Thus in this case 
the dimension of the space Bi will be two. Hence there is no hope of getting the inf-sup 
condition for this choice of the bubble function. 

4.1.2 The first choice of bubble functions 

In the next step, we consider the bubble function 


b K = &A<p K xy{l -x)(l-y), 
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where tp K is the standard bilinear basis function corresponding to the lower-left corner of 
the square K. Since ipx = (1 — x)(10?/), the bubble function bn on the reference square I\ 
can be defined as 

b K = 64(1 - z)(l - y)xy{ 1 - x)(l - y). 

Defined in this way the bubble function bx does not depend on the local numbering of the 
vertices of K. In this case, the matrix D has rank 8, and is computed as 
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Remark 2 We have used the gradient of the bilinear function ipx to construct a vector 
bubble function associated with the element K in m- Since the construction of the bubble 
function using the gradient of tpK cannot be done on a reference element, this new bubble 
function is computationally much easier. 

4.1.3 The second choice of bubble functions 

It is interesting to see if we can multiply the bubble function by a linear function and obtain 
the stability. For this purpose we can choose a bubble function on the unit square (0, l) 2 as 

bK = (a + bx + cy)xy{ 1 — a;)(l — y), abc ^ 0. 

For simplicity we choose 

b K = 8(1 + x + y)xy( 1 - x)(l - y). 

We note that the factor 8 is used to force the value of the bubble function at the centroid 
of the square to be 1. The resulting matrix D has also rank 8 in this case, and hence the 
dimension of the space Bi is one. Moreover, the matrix D is computed as 
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D = 












Remark 3 The proof of stability is presented for the two-dimensional case. However, this 
can be extended to the three-dimensional case without a major change. 

Remark 4 It is interesting to see if we can use a quadratic function symmetric about the 
centroid of the element to multiply the standard bubble function. To check this we use a 
bubble function on the reference square K = (0, l) 2 defined as 

b K = xy (^x 2 + y 2 - x - y + y ^ (1 - x)(l - y), 

and compute the matrix D. In this case, the rank of the matrix D is just 1, and hence the 
dimension of the space B.i is 2. 


An immediate consequence of the above discussion is the well-posedness of the discrete 


problem (3.4). From the theory of saddle point problem, see, e.g., 0, we have the following 
theorem. 


Theorem 5 The discrete problem (3.4) has exactly one solution ( Uh,Ph ) £ F^x S^, which 


is uniformly stable with respect to the data f, and there exists a constant C independent of 
the mesh-size h such that 


||w/i||i,fi + |b/Jo,n < C||/||o,n- 


The convergence theory is provided by an abstract result about the approximation of saddle 
point problems, see [5]. 


Theorem 6 Assume that ( u,p ) and (uh,Ph) be the solutions of problems (2.2) and (3.4), 
respectively. Then, we have the following error estimate: 


u h \\i,n + lb ^ Ph ||o,o < C 


inf 
v h eV h 


\\u-v h \\ lyn + 


inf lb-9/Jo,n 


(4.2) 


5 Numerical Results 

In this section we present two numerical experiments to verify the optimal a priori error 
estimate and some numerical experiments to verify the inf-sup condition for the proposed 
finite element scheme. For both examples we consider a simple unit square fl = (0, l) 2 . 


5.1 Verify a priori error estimate 

For both examples we consider a uniform initial triangulation consisting of four squares. 

First example. For the first example we choose the exact solution u = (ui,u 2 ) as 
Mi = -2 x 2 y (2 y - 1) (x - l) 2 {y - 1), u 2 = 2 x y 2 (2 x - 1) (x - 1) (y - l) 2 . 

We use the kinematic viscosity v = 1. The exact solution for the pressure is chosen as 

P = x{ 1 - x){l - 2 y), 

so that p G Lg(fl). The exact solution u satisfies the homogeneous Dirichlet boundary con¬ 
dition on dkl, and the right hand side function / is computed by using the exact solution u 
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and the pressure p. We have presented the errors in the velocity and the pressure approx¬ 
imation using the H 1 -norm and the L 2 - norm, respectively in Table [l] for the first choice 
of the bubble function, and in Table [2] for the second choice of the bubble function. We 
note that the standard choice of the bubble function leads to a singular matrix. From the 
presented tables we can see the optimal convergence of the velocity approximation in the 
H 1 and L 2 -norms, and a super-convergence result for the pressure in the L 2 -norm. As we 
expect a convergence rate of order 1 for the pressure approximation in the L 2 -norm but get 
a better approximation of order 1.5, this is a super-convergence. This better convergence is 
due to the fact that we have used the standard continuous bilinear finite element space for 
the pressure approximation. We can also observe that all errors are smaller for the second 
choice of bubble functions. 


Table 1: Discretization errors for the velocity and pressure, Example 1 (First choice) 


level l 

# elem. 

lb - uji.n 

q 

o' 

53 

1 

53 

lb-pjo.fi 

1 

16 

3.23129e-02 


3.03116e-03 


1.76150e-02 


2 

64 

1.58286e-02 

1.03 

8.24246e-04 

1.88 

7.00356e-03 

1.33 

3 

256 

7.79938e-03 

1.02 

2.06421e-04 

2.00 

2.50753e-03 

1.48 

4 

1024 

3.87699e-03 

1.01 

5.12144e-05 

2.01 

8.78516e-04 

1.51 

5 

4096 

1.93346e-03 

1.00 

1.27289e-05 

2.01 

3.08875e-04 

1.51 

6 

16384 

9.65545e-04 

1.00 

3.17131e-06 

2.00 

1.08856e-04 

1.50 


Table 2: Discretization errors for the velocity and pressure, Example 1 (Second choice) 


level l 

# elem. 

||w - wji.fi 

q 

0 

s 

1 

S3 

IIP -Ph||o,fi 

1 

16 

3.16876e-02 


2.89325e-03 


1.17765e-02 


2 

64 

1.56503e-02 

1.02 

7.90369e-04 

1.87 

4.31789e-03 

1.45 

3 

256 

7.75922e-03 

1.01 

1.99983e-04 

1.98 

1.44890e-03 

1.58 

4 

1024 

3.86716e-03 

1.00 

4.99365e-05 

2.00 

4.93948e-04 

1.55 

5 

4096 

1.93102e-03 

1.00 

1.24544e-05 

2.00 

1.71287e-04 

1.53 

6 

16384 

9.64934e-04 

1.00 

3.10849e-06 

2.00 

5.99594e-05 

1.51 


Second example. For the second example we consider an exact solution given in [3], 
where the exact solution for the velocity u = (ui,U 2 ) is given by 

ui = x + x — 2 xy + x — dxy + x y, 112 = — y — 2 xy + y — dx y + y — xy , 

and the exact solution for the pressure is given by 

32 4 

p = xy + x + y + xy - 

We use the kinematic viscosity v = 1 and the exact solution to compute the right-hand side 
function /. As in the first example we compute the errors in the velocity and the pressure 
approximation using the H 1 - norm and the L 2 - norm, respectively. The numerical results 
are tabulated in Table [3] and [4] for the two choices of bubble functions, respectively. As in 
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the first example, we can see the optimal convergence rates for the velocity approximation 
in H 1 and A 2 -norms, and a better convergence rate for the pressure in L-norm. We also 
observe that all errors are smaller for the second choice of bubble functions although the 
difference is quite small in this example. 


Table 3: Discretization errors for the velocity and pressure, Example 2 (First choice) 


level l 

# elem. 

||« - w/i||i,n 

||« - UfcHo.n 

IIP-Philo,n 

1 

16 

6.96126e-01 


3.33821e-02 


2.25132e+00 


2 

64 

3.39100e-01 

1.04 

8.37772e-03 

1.99 

5.58680e-01 

2.01 

3 

256 

1.66684e-01 

1.02 

2.09556e-03 

2.00 

1.59539e-01 

1.81 

4 

1024 

8.26546e-02 

1.01 

5.24458e-04 

2.00 

4.49273e-02 

1.83 

5 

4096 

4.11633e-02 

1.01 

1.31193e-04 

2.00 

1.28191e-02 

1.81 

6 

16384 

2.05425e-02 

1.00 

3.28081e-05 

2.00 

3.80370e-03 

1.75 


Table 4: Discretization errors for the velocity and pressure, Example 2 (Second choice) 


level l 

# elem. 

||m - % l,f! 

||w - UfcHo.n 

IIP-Philo,n 

1 

16 

6.96024e-01 


3.23184e-02 


5.93926e+00 


2 

64 

3.35337e-01 

1.05 

7.82819e-03 

2.05 

4.04732e-01 

3.88 

3 

256 

1.65795e-01 

1.02 

1.97572e-03 

1.99 

6.07983e-02 

2.73 

4 

1024 

8.24467e-02 

1.01 

4.97135e-04 

1.99 

1.78268e-02 

1.77 

5 

4096 

4.11137e-02 

1.00 

1.24714e-04 

2.00 

5.88206e-03 

1.60 

6 

16384 

2.05304e-02 

1.00 

3.12328e-05 

2.00 

1.98964e-03 

1.56 


6 Conclusion 

In this contribution we present a finite element method for Stokes equations using continuous 
bilinear finite elements enriched with bubble functions for the velocity approximation and 
continuous bilinear finite elements for the pressure. In contrast to an earlier contribution we 
show that a single vector bubble function per element is enough to guarantee the stability of 
the discrete linear system. The numerical results also demonstrate the optimal convergence 
rates for the velocity and pressure approximation. 
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